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Abstract 

One intriguing beyond-the-Standard-Model particle is the QCD axion, which could simulta¬ 
neously provide a solution to the Strong CP problem and account for some, if not all, of the 
dark matter density in the universe. This particle is a pseudo-Nambu-Goldstone boson of the 
conjectured PecceRQuinn (PQ) symmetry of the Standard Model. Its mass and interactions are 
suppressed by a heavy symmetry breaking scale, fa, whose value is roughly greater than 10® GeV 
(or, conversely, the axion mass, ma, is roughly less than 10^ yeV). The density of axions in the 
universe, which cannot exceed the relic dark matter density and is a quantity of great interest in 
axion experiments like ADMX, is a result of the early-universe interplay between cosmological evo¬ 
lution and the axion mass as a function of temperature. The latter quantity is proportional to the 
second derivative of the temperature-dependent QGD free energy with respect to the CP-violating 
phase, 9. However, this quantity is generically non-perturbative and previous calculations have 
only employed instanton models at the high temperatures of interest (roughly 1 GeV). In this and 
future works, we aim to calculate the temperature-dependent axion mass at small 9 from first- 
principle lattice calculations, with controlled statistical and systematic errors. Once calculated, 
this temperature-dependent axion mass is input for the classical evolution equations of the axion 
density of the universe, which is required to be less than or equal to the dark matter density. 
Due to a variety of lattice systematic effects at the very high temperatures required, we perform a 
calculation of the leading small-0 cumulant of the theta vacua on large volume lattices for SU(3) 
Yang-Mills with high statistics as a first proof of concept, before attempting a full QGD calculation 
in the future. From these pure glue results, the misalignment mechanism yields the axion mass 
bound ma > (14.6 ± 0.1) yeN when PQ-breaking occurs after inflation. 

PACS numbers: 12.38.Gc, 14.80.Va 
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I. INTRODUCTION 


Despite overwhelming experimental and theoretical evidence that Quantum Chromody¬ 
namics (QCD) is the underlying theory of interactions between quarks and gluons, one 
puzzle has eluded explanation for over 30 years: the Strong CP problem [TH3] . QCD allows 
for a 0{1) combined charge-conjugation and parity (CP) violating phase 9 and, yet, this 
phase is found experimentally to be consistent with zero to within one part in ten billion |1] . 

Other than attributing this result to an appreciable hne-tuning, there have been three 
proposed explanations explored at length: the possibility of a zero up-quark mass [5] , spon¬ 
taneous CP breaking [6HT2], and an additional “hidden” chiral symmetry. After taking 
into account the up-to-date Standard Model flavor-changing constraints [13] and lattice 
QCD calculations [T3|, the most feasible proposal is the last. The additional chiral symme¬ 
try proposed, also referred to as Peccei-Quinn (PQ) symmetry [T51 [12], is such that upon 
spontaneous symmetry breaking, the effective potential naturally gives a zero CP-violating 
phase^. As with any spontaneous breaking of a continuous symmetry, one would expect a 
resulting (pseudo-)Nambu-Goldstone particle to be present in the universe; this particle is 
known as the axion^ HaUHlEIlES]. 

The original proposals focused on axions masses at or below the electroweak scale [I5l |T6] , 
but collider [23l |23] and astrophysical constraints [251 - 127] now require the axion to have a 
mass below 10^ /reV. A suitable class of axion models, named “invisible axions”, could 
allow for light axions whose interactions are suppressed by very high energy scales [2HH3I]- 
In these models the large density of light axions could potentially account for some, if not all, 
of the dark matter density in the universe [32H31]- For this reason, dedicated experiments 
such as ADMX [351 - 137] search for the axion coupling to two photons directly. With the 
next generation of axion experiments underway 1371 along with increased constraints on 
inflation [38], a great deal is expected to be learned about potential axion interactions in 
the next few years. 

There has been a wealth of research on axion cosmology and the related axion energy 
density over the past 30 years [32H311 EHllSS] • One particular constraint, the “overclosure 
bound”, requires that the axion density today not be greater than the total observed present- 
day dark matter abundance. This bound requires accounting for the evolution of the axion 
mass through cosmological history and its consequences for axion production. Depending 
on when PQ breaking occurs relative to inflation, the overclosure bound gives different 
relations. If PQ breaking occurs during or before inflation, the bound is on a relation of 
two variables, the initial value for the CP-violating phase inside our cosmic horizon and 
the axion decay constant, fa- If PQ breaking occurs after inflation, the bound is solely on 
fa- After inflation, there are three stages to this early universe, high-temperature evolution: 
the post-inflation evolution to a time when the axion mass is comparable to the inverse 
horizon of the universe (i.e. the Hubble constant), the subsequent evolution in the chirally 
symmetric phase, and the time period before and after the chiral symmetry breaking phase 
transition of QCD. In the hrst stage, the axion mass depends greatly on the QCD free energy 
which is a function of the CP-violating phase, 9, and temperature, T. To date, the only 
methods employed for estimating this QCD free energy dependence as related to the axion 


^ Even this explanation for a small 6 has limitations, since PQ-violating Plank-scale operators up to dimen¬ 
sion 10 can reintroduce the Strong CP problem [ndH]. 

^ For some general reviews on the subject of axions, see Refs. [MHO]. 
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mass are the dilute instanton gas model (DIGM) and the interacting instanton liquid 
model (IILM) [58]. The DIGM is expected to be valid only for high temperatures, while the 
IILM models strong interactions between instantons around the QGD phase transition. In 
the context of axions, both of these scenarios have been explored in detail [STl I5U] . While 
the desired characteristic temperatures are indeed large (~ 1 GeV), one open question is to 
what extent the DIGM and IILM are valid (i.e. to what extent non-perturbative physics 
plays a role at these temperatures and is modeled properly) and how accurate the overall 
scale is (the free energy scales as the conhnement scale, Aqcd, to a large power). Another 
related issue is that neither the DIGM nor the IILM yield controlled uncertainties from hrst 
principles. The only known approach that can address these two quandaries^ quantitatively 
is lattice gauge theory, which is what we apply in this work. 

There have been ample lattice studies of QGD vacuum properties at high temperature, 
highlighted by recent work verifying the QGD critical temperature {Tc) with 2+1 quarks at 
physical pion masses in a chiral discretization [60], and controlled continuum limits for the 
QGD equation of state [611162] with temperatures as large as T ~ 600 MeV. Unfortunately, 
simulating higher temperatures becomes computationally more expensive due to unphysical 
systematic hnite-volume effects (for a hxed number of lattice sites, the physical volume 
decreases as the temperature increases). However, simulations of pure Yang-Mills with 
three colors is appreciably cheaper and high-temperature 6 = 0 studies spanning T ~ 
(5 —1000)Tc suggest that non-perturbative effects can still be appreciable at T ~ 1 GeV [63] • 
Appreciably more difficulty ensues when exploring 6* fy 0 quantities, as the topological term 
in the Lagrangian has an associated sign problem that renders standard lattice Monte- 
Garlo techniques intractable^. However, techniques do exist for extracting quantities at 
small 6, such as the leading 6'^ dependence of the free energy which is proportional to the 
topological susceptibility®. While topological fluctuations are a topic of active focus for 
lattice calculations at zero-temperature [6MQ], there have only been a handful of studies at 
finite temperature; first for Yang-Mills and more recently in full QGD using chiral 

discretizations eg. 

In this work, we aim to extend and improve the finite-temperature results in Ref. [72] with 
the express purpose of comparing with the DIGM and IILM and, within the capability of this 
calculation, quote a first-principles bound for the axion mass with controlled uncertainties. 
It is important to note that this work is only an initial step towards the full QGD problem 
and it contains two primary inadequacies due to limited resources and lattice technology. 

The first inadequacy is that our calculations are performed for a 3-color pure Yang-Mills 
theory without the dynamical fermions of QGD. Ultimately, this may not prove to be too 
far from the physical result, as effects of fermion loops are expected to be suppressed at 
high temperatures [76]. For topological quantities, however, this may not be the case m 
and full QGD calculations should be pursued in the future. The reasoning for studying a 
purely gluonic theory at this stage is two-fold. First, Yang-Mills theories are over ten times 
cheaper to simulate than full QGD due to the ability to employ the heatbath algorithm 
for entirely bosonic theories. This extra gain in computational efficiency will allow us to 
go to higher temperatures and larger volumes, and to directly estimate systematic effects 
for topological quantities. The second advantage is that heatbath Monte-Garlo calculations 


^ For these reasons, the use of lattice QCD was suggested in Ref. [51] . 

^ For progress in solving sign problem for actions with non-zero 0, see Ref. |64L 165] . 
® For a review on lattice calculations of topological susceptibilities, see Ref. [55] . 


3 



have appreciably shorter autocorrelations in extracting topological quantities than the hybrid 
Monte-Carlo algorithm used for dynamical fermions. This will allow us to extract the high 
level of independent statistics required to address the problem at hand. 

The second inadequacy in our present calculation is that reaching the large values of 
T/T, 5 with controlled volume systematics is still beyond our current computational re¬ 
source limitations®. Summarizing our calculation, we aim to test the validity and overall 
scale of the DIGM/IILM for a SU(3) Yang-Mills theory for T ~ 372 — 710 MeV (or, alter¬ 
natively, T/Tc ~ 1.31 — 2.50, remembering that Tc for SU(3) Yang-Mills is ~ 284 MeV [78]) 
with high statistics and controlled volume and discretization uncertainties. From these re¬ 
sults, we extrapolate our results to extract the characteristic temperature when the axion 
mass and Hubble constant are comparable and proceed to evaluate what this result would 
imply for the axion mass in the present day universe with propagated uncertainties. 


The next section is a brief review of the cosmological evolution of axions, while Sec. Ill 


describes the free energy of QCD in the presence of a 6* term and a hnite temperature. 


The paper then continues with a detailed description of the lattice simulations in Sec. IV 


lattice results in Sec. |V] and lattice error budget in Sec. |V^ Based on hrst-principles lattice 
results we derive a bound on the axion mass from the aforementioned overclosure argument 
in Sec. VII and Sec. VIII before our concluding remarks. 


II. BRIEF REVIEW OF AXION COSMOLOGY 


The theory of quarks and gluons. Quantum Chromodynamics (QCD), supports a term in 
its action iSqcd which violates the combined charge-conjugation and parity symmetry (CP), 


‘5qcd 3 dQ 


Q= d'^x 


327r2 


tiFF 


( 1 ) 


where Q is the topological charge, F (F) is the (dual) gauge held strength tensor, and g 
the QCD gauge coupling constant. The Strong CP problem is the observation that the 
parameter 6 could in principle take any value between the maximally-CP-violating values 
of —71 and TT, but is measured to be consistent with the CP-conserving value of zero to 
one part in ten billion |3]. An elegant explanation of the small value of 6 was proposed by 
Peccei and Quinn: it is possible to promote 0 to a dynamical variable in such a way that 
it is controlled by a (pseudo-)Nambu-Goldstone boson called the axion. The Peccei-Quinn 
(PQ) symmetry makes the axion naturally light and creates a potential for 6 which favors 
small 6. The QCD-|-axion Lagrangian is given by 


9 


mg)q + 


327r2 




T^a T^a{iu 


( 2 ) 


where fa is proportional to the vacuum expectation value that breaks PQ symmetry. This 
quantity is the one free parameter of standard QCD axion theories, and it is this scale which 
ultimately determines the axion mass and its couplings to two photons that experiments are 
actively searching for [3^1 - [37] . 


® Future explorations of extracting topological quantities on anisotropic lattices may alleviate some of the 
computational burden on volume as long as discretization effects stay minimal. 
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While fa is the primary parameter for both the axion mass and axion energy density 
given the QCD Lagrangian, non-perturbative QCD effects lead to non-trivial cosmological 
conseqnences, especially for the axion nnmber density. At low temperatnres, after chiral 
symmetry has been broken, the axion mass and conpling obeys a relatively simple rela¬ 


tion 


rriafa 


y/murud 
rriu + rud 




(3) 


where m.,t and f-^ are the pion mass and pion decay constant, respectively. The axion mass 
at temperatnres in the chirally-symmetric phase of QCD will be the primary focns of this 
work. 

One signihcant constraint on axions is that their present energy density does not exceed 
the dark matter density of the nniverse (the “overclosnre” bonnd). To determine the rel¬ 
evant energy density today, one must explore the interplay between the axion mass and 
cosmological evolution by solving the equations of motion. However, the initial question in 
this evolution is to place the moment PQ symmetry breaks along the cosmic timeline. There 
are two options: 


1. PQ symmetry breaks before or during inflation. 


2. PQ symmetry breaks after inflation. 

While the hnal energy density comes down to one (or effectively two) free parameters, 
there are many non-trivial calculations and subtleties that need to be understood from both 
cosmology and non-perturbative held theory. It is useful to give a brief summary on how 
the energy density arises in each of these PQ-breaking scenarios^: 

1. For PQ symmetry breaking before or during inhation, the held 9 = a/fa would be 
homogenized over large distances. In other words, the causally disconnected regions of 
spacetime that have diherent initial values of the held (angle) before inhation stretch 
to cosmic scales, so that our universe has a uniform initial 9. Any excited axion 
mode or topological defect (e.g. strings) will be diluted away, leaving only the zero¬ 
mode contributions to the energy density. The resulting axion density arises from 
this “misalignment mechanism” and has large dependence on this initial theta angle 
(often referred to as the misalignment angle, 9i, which can take any value between 
0 and 27r). Ultimately, the hnal energy density will be proportional to roughly 
which is ehectively a free parameter along with fa- Thus, the overclosnre bound can 
only bound fa as a function of 9i. Also, due to large huctuations during inhation, 
isocurvature density perturbation bounds from CMB observations can put important 
constraints on 9i as a function of the Hubble scale during inhation [561 ESHU] . 

2. For PQ symmetry breaking after inhation, the misalignment angle 9i is (ehectively) 
averaged over, since all values G [0, 27r] are equally probable in small regions of our 
universe. Moreover, there are several additional ehects that must be accounted for 
in addition to the misalignment mechanism in this case. First, nonzero-momentum 
modes can contribute and their ehect is discussed in Refs. ga Eg. Second, when 
PQ symmetry is broken, global topological defects called axionic strings will form 


^ For detailed reviews, axion cosmology, see Ref. [HI EH [56] 
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and decay to axions, introducing an efficient process for energy loss and effectively 
increasing the total axion number density. It is currently debated as to how much 
these decays will increase the total axion number, with some claiming between one and 
two orders of magnitude more than the misalignment mechanism [101 021 03 EDI 05] , 
while other claim it is on the same order 53 03 ESj. The third effect is the axion 
strings can be connected to one or more domain walls EH which correspond to different 
minima of the axion potential. The decay of these domain walls also leads to additional 
axions, but is believed to be subdominant to string decay [13EHI ES] . Overall, the 
energy density can be calculated in terms of one free parameter, /„ and, as a result, 
an overclosure bound on fa can be derived. 


Our goal in this and future work is to improve the non-perturbative QCD input that goes 
into the axion density calculations. In particular, we aim to provide a controlled calculation 
of the temperature-dependent axion mass from hrst-principle lattice calculations®. The axion 
mass in the chirally symmetric phase is given by 


ml{T)fa = 


d^F{e,T) 


de^ 


^ x{T), 


e=o 


(4) 


where F{6, T) is the QCD free energy as a function of CP-violating phase and temperature, 
while X is the topological susceptibility. In each of the energy density scenarios discussed 
above, the temperature-dependent axion mass ma{T) plays a role, particularly when this 
quantity is comparable to the Hubble scale in the early universe evolution. We will restrict 
our discussion primarily to the analytic evaluations of the misalignment mechanism, as the 
lattice calculations performed in this work could be appreciably different from QCD due the 
computational and algorithmic limitations discussed in subsequent sections. While we will 
not discuss it in detail, we will also point out how the axion mass enters into the relevant 
calculations for axion strings and domain walls. Once the full QCD results are at a mature 
stage, the lattice result for rriaiT) should be used as input for the numerical solutions to the 
classical equations of motions and cosmology simulations. 


A. Misalignment Mechanism 

To keep the presentation self-consistent, we summarize the description of the misalign¬ 
ment mechanism of Ref. [32] and subsequently summarized in Ref. [521 IHTl l56] . We start 
from from the Robertson-Walker metric of the universe, 

— ds^ = dt^ — R{tYdx ■ dx , (5) 

and write down the axion equation of motion: 

+ 3Hdt - aix) + ml{t)fa sin = 0 , (6) 


® Other non-perturbative QCD aspects include the number of theta-vacua/domain walls for 9 changing by 
27r [57] and whether or not axion number density is constant throughout the QCD phase transition [551155] . 
We do not discuss these further in this work. 
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where H = R/R is the Hubble constant and the second term is given by the derivative of 
the effective potential for the axion held, 14 , given by. 




1 - cos I — 
Ja 


(7) 


Along with the axion equation of motion, the Hubble constant evolution is given by the 
Einstein equation 


= 


1 


TT' 


1 f da\ 


3M2 


+ +rnl(t) 


30 


dt J 


cos 


fa 


( 8 ) 


where Mpi is the reduced Planck mass ^/hc/SirG. For given ma{t), or equivalently ma{T), 
one can solve Eq. @ and Eq. ([^ numerically to arrive at the axion energy density: 


Pa 


1 / da\ 




(9) 


At this point, it is useful to discuss the qualitative consequences of these solutions (see 
Fig. 4 in Ref. [51j for an illuminating plot). When considering early times, one notes that 
the Hubble constant is much larger than the axion mass (3if 3> rria)-, meaning that the 
axion wavelength is larger than the Hubble length. Thus, the axion does not feel a potential 
and it is effectively massless. Moreover, the axion number density is zero and the axion 
held has a constant value, di (the misalignment angle). As time increases (and temperature 
decreases), the axion mass increases while the Hubble constant decreases. The solutions 
change drastically when the axion mass is of the same order as the Hubble constant (3if ~ 
ttIq), at which point the axion mass “turns on” and the axion held rolls down the potential— 
during this period the axion number density jumps to a nonzero value. From this point in 
time onward, both the axion held and axion number density have decaying oscillations about 
their eventual hnal values that we see today (the decay becomes adiabatic when 3iif ^ rn^). 

Analytic progress can be made by making a few key observations®. First, let us quantify 
the characteristic temperature and corresponding time, Ti and ti, respectively, when the 
Hubble constant is comparable to the axion mass 

m^Ti) = 3hf(Ti) = ^ 

When T >Ti, Eq. (|^ reduces to 

7r^y^(T)T^ 

y/^Mpi 



( 11 ) 


( 10 ) 


where g^{T) is the ehective number of relativistic degrees of freedom at a given temperature. 
In this work, we use the parameterization of g^{T) in Appendix A of Ref. [5l]- We will 
compare our lattice niaiT) to H{T) using Eq. (10) to extract Ti. The second observation is 


9 


For simplicity, we look at just the zero-momentum mode, dropping the from Eq. §. For a complete 
calculation with the non-zero modes, see Ref. [45115^ . 
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that the decaying oscillations when 3H nia are adiabatic and the eqnation of motion can 
be recast as 


1 dpa ^ ,^.dma 


1 - cos I — 
Ja 


dt 


/2 dt dt 

When rUa 3> H and S> drUa/dt, this relation leads to the adiabatic invariant 

PaR^ 


rUn 


= constant , 


( 12 ) 


(13) 


and snbseqnently [M] . 


PaiT^) = pa{Ti) 


rUaiT^) V paiTi) m„(T^) 9.iT^)T^ 


ma{Ti) \Rm 


7 m^iT,) g,{T,)Tr 


(14) 


where s is the entropy density, ~ 2.73 K is the temperatnre of the cosmic microwave 
backgronnd today and 7 is the ratio of the entropy density today to the entropy density at 
ti. From this relation, along with the fact that Pa{Ti) ~ and Eq. Eq. (10), 

and Eq. 0 . we arrive at the commonly used relation 


PaiT^) = 


^/90g^{Ti) mu + md Mpi 


fa 

Ti 


FijOi 

27 


eh 


(15) 


where Fi{6i) accounts for the anharmonic corrections to Eq. (13). If it is also assumed 


that the expansion of the universe is adiabatic for temperatures below Ti, the quantity 
Fi( 6 'i )/27 1 for 6^1 < 2 [53]. Using lattice QCD calculations, Ti can be extracted as a 

function of fa and a conhnement-scale observable (such as the QCD deconhnement tem¬ 
perature Tf), and as a result, a bound on the free parameters fa and 9i can be derived 
to ensure that the axion density is below that of the dark matter density. If we are only 
exploring axion theories where PQ-breaking occurs after inflation, the misalignment angles 


are averaged over, and the 9i dependence in Eq. (15) is replaced with the average value, 
ef -t {9j) = n/y/3. 


B. Axion Strings 

In the case where PQ breaking occurs before or during inflation, any resulting topological 
defects, such as axion strings or domain walls, are diluted away and only the misalignment 
mechanism contributes to the axion density with the hxed 9i of our observed universe. 
However, if PQ breaking occurs after inflation, these topological defects can decay into 
axions and signihcantly alter the number density of axions in the universe today. Scenarios 
with two or more decaying domain walls are highly constrained from neutron EDM bounds 
on CP violation and would require an appreciable hne-tuning |5B] . 

The current understanding is that the quantity of axions produced by domain wall decay 
is below the number of axions produced from string decay, a number which is at least 
comparable to the density produced from the misalignment mechanism and could even be 
appreciably larger as debated in the literature [4011421147[ I49| l5n[ l85| l 86 ] . In the case of string 
decay, this process occurs for string frequencies between the Hubble scale and the axion mass 
{H < u < ma), and thus begins to play a role at temperature Ti. Similarly, domain walls 



















also have a dependence on rua, which plays a role in the axion spectrum. While it is not 
generally expected to be the primary source of uncertainty in these large scale cosmological 
simulations [Ml ES] , accurate input for the axion mass as a function of temperature would 
nonetheless be useful. 


III. QCD FREE ENERGY AS A FUNCTION OF T AND 9 


The free energy of QCD, F, as a function of temperature and theta is given in terms of 
the path integral 


Zqcd{0,T) = I [dA][d'ijj][d'ijj] exp ( -T '^d^x Cqcd{6) ) = exp[-UF(6',T)], (16) 


where 

Cqcd(0) = Cqcd + (17) 

The free energy is periodic in 9 and thus we can restrict to the range —tt < 6 < n and 
parameterize the free energy in terms of a sum of cosine functions, 

F{9,T) = Y,Cn{T)cos{n9). (18) 


Since QCD is a non-perturbative problem for generic T and 6, the values of Cn{T) are not 
readily accessible; even in hrst-principle lattice QCD calculations, non-zero 9 introduces a 
computationally intractable sign problem. However, at high enough temperatures, QCD 
interactions are perturbative and their parametric dependence on T and 9 should be given 
by the dilute instanton gas model (DIGM) [57] • In this model, where the QCD background 
is approximated by non-interacting instantons, only Ci contributes to the free energy (all 
n > 1 coefficients are zero). The DIGM also predicts the value of the coefficient (whose 
mass dimension is 4) as a function of Aqqd, 

Fdigm{9,T) = -C{T)cos{9), (19) 


where C (T) can be well-approximated by [53l ES] 


G(T)^ 


CA^ 

(T/A)"’ 


( 20 ) 


where G is a dimensionless constant in temperature and the fermion masses and n is a 
power that is a function of number of flavors and colors^°. The A in this equation essentially 
represents the scale setting of QGD, which will be key topic of discussion later in this paper. 
For three flavors, the latest value found is G ~ 1.274 x 10“^^ for A = 440 MeV [53 l54] . 
However, this value can vary by almost an order of magnitude if A is varied by 15%. 

The interacting instanton liquid model (HLM) is a more sophisticated model that ac¬ 
counts for instanton-instanton interactions and was applied to the problem at hand in 


The actual expression derived from DIGM contains a more involved integral over the logarithmic running 
of the coupling. For a complete expression see Appendix A of Ref. [M] 
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Ref. [58]. To solve for the free energy in this model, grand canonical Monte-Carlo simu¬ 
lations of the partition function are required. The results can be fit to the form 

FuLM{e,T) = -D{T)cos{e), ( 21 ) 


with 


D{T) ~ 


grfoA4 

(T/A)-'^i 


exp 



+ da 



( 22 ) 


where the values for di are given in Ref. [5l| as a piecewise function of temperatures at the 
flavor mass thresholds (coincidentally, values for do and di do not differ from the analogous 
terms calculated in Ref. [5S|)- The authors point out that the largest uncertainties arise 
from the scale setting of using A = 400 MeV and emphasize at various stages that the IILM 
should be compared with lattice QCD results. 

Since lattice QCD calculations are numerical and ultimately yield dimensionless numbers, 
lattice scale setting is of vital importance to relate to lattice results to reality. The way this 
is typically done is to choose an observable (preferably one that can be calculated to high 
precision with little to no unphysical lattice artifacts from volume or hnite lattice spacing) 
and match it onto some measured/derived quantity from experiment. This is often done with 
matching onto heavy-quark potentials [HnHHI] , but this can also be done for spectroscopy of 
conhnement scale masses, such as the omega baryon [22] • For our pure-glue calculation, the 
most reliable scale setting is to use the string tension, a, which can be used to dehne the 
critical temperature Tc in physical units. Once this is done, all scales can be expressed in 
units of the critical temperature Tc which makes it natural to fit to the forms 


CT} 

iT/K) 





+ da 



(23) 


It should be noted that when using the string tension (or any other heavy-quark scale) for 
setting the scale, the critical temperature in Yang-Mills [TS] is roughly twice that of full 
QCD [Sni ESI El] • This is to be expected, as the dynamical quarks play a non-trivial role in 
the phase transition. 


IV. LATTICE SIMULATIONS AND SCALE SETTING 

We perform the lattice calculation of the free energy and its expansion in the 6 angle 
using the Wilson plaquette action to discretize the SU(3) Yang-Mills continuum theory: 


J]] i 1 - -Re tr[t;p) 


(24) 


where Up is the ordered product of gauge links along the plaquette and (3 is the lattice 
gauge coupling. This lattice action is well known and, in particular, a great deal of data 
exists in the literature regarding the scale setting procedure of the corresponding lattice 
system. We will make use of this information when setting the temperature scale of our 
lattice simulations. The lattice action in Eq. (24) contains only one free parameter, f3, 


which sets the lattice spacing through dimensional transmutation. Once the lattice spacing 
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is fixed, one can give it physical nnits by measuring a physical observable and relating it to 
an experimental value. 

We are interested in spanning a large interval of temperatures with our simulations. In 
order to accomplish that we £x the number of lattice sites in the temporal direction to 
6 and set the physical temperature 


T(/3,W) = 


a{/3)K 


(25) 


by tuning a value for the lattice coupling constant f3. In Eq. (25) we explicitly show the 


dependence of the lattice spacing a on jS. Working at a fixed W, the lattice spacing alone 
determines the temperature. This approach is widely used in the literature [63]. To check 
for possible lattice discretization effects in our results and establish a continuum limit, 
we simulate the same temperature at a smaller lattice spacing by simply increasing the 
temporal lattice sites (we use = 8). Once the lattice spacing is fixed, the physical spatial 
volume is determined by the number of lattice points W in the three spacial dimensions: 

= {aNcr)^. One complication of our approach is that the lattice spacing gets smaller at 
higher temperatures, implying smaller spacial volumes at hxed aspect ratio N„/Nt-. As a 
consequence it is of paramount importance to check for possible hnite-volume effects at high 
temperatures. 


Following Eq. (25), setting the physical temperature scale has been traded for deducing 


the lattice spacing as a function of the bare lattice parameters. In order to do so, one can 
follow a variety of approaches that usually involve measuring a physical quantity with the 
dimensions of a mass or a length in a zero-temperature setup. We adopt the strategy of 
using the string tension a, which is well understood for a Yang-Mills theory: 


T _ a^{(3c)N^c 

Te a^{(3)Nr ^ ’ 

where (3c gives the critical temperature Tc in a box of thermal extent Nt-c from the temporal 
Polyakov loop susceptibility. In particular we choose the most continuum-like point at which 
the thermal phase transition has been studied for SU(3) with Wilson plaquette action [nS] 
{f3c = 6.338, Nrc = 12) and we interpolate numerical results for a^/a as a function of (3 from 
Ref. [95]. The interpolation uses the same function described in Ref. I7HI and later rehned 
in Ref. [96]. In Fig. [^we show T/Tc{(3) for two temporal extents W = 6, 8. 

We can compare this approach with scale setting methods that use a different physical 
quantity, like the Sommer radius vq, and quantify a systematic error for our procedure. 
We have verihed that setting the scale via both the static quark potential quantity tq, as 
described in Ref. ra. and the continuum-extrapolated ratio of the critical temperature and 
the string tension Td\f(j as described in Ref. [HS] give T/Tc values which agree with our 
method up to 1% - 1.5% corrections in the temperature range explored in the paper. Across 
the rest of the paper we report either results in units of the lattice spacing or results in 
units of the critical temperature. This is to avoid giving physical units to a or Tc; for a 
SU(3) Yang-Mills theory like the one we simulate, the systematic error associated to giving 
physical units to the lattice quantities is of the 4% - 7% level [69] . 

Our gauge configurations are generated using the Chroma software system [98] with GPU 
acceleration supplied by QDP-JIT/PTX [HH]- The update algorithm is a standard Cabibbo- 
Marinari heatbath HM] for SU(3), alternated with 8 steps of standard over-relaxation to 
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FIG. 1: Temperature as a function of /3, given /3c = 6.338 and Nr = 12. 


reduce correlations among subsequent configurations. For each temperature T/Tc we sim¬ 
ulate two complementary streams of configurations: one starts from a random gauge con¬ 
figuration, while the other is initialized to a conhguration where each link is a unit matrix. 
On both streams we discard the hrst 10000 updates and monitor basic local observables like 
the spatial and temporal plaquettes and the Polyakov loops in the four dimensions. For all 
temperatures the two streams thermalize to a common value of plaquette and Polyakov loop; 
this allows us to effectively double our statistics. More details on the topological charge and 
topological susceptibility analysis are given in the next section. Our full set of ensembles 
and measurements is shown in Table [H 

V. TOPOLOGICAL SUSCEPTIBILITY ON THE LATTICE 

We measure the topological charge on the lattice using the bosonic held-theoretical deh- 
nition that is built upon discretizing the continuum formula on the right of Eq. Q; 

X 

where U^y{x) is the /xz/-plane plaquette at the lattice point x. The above dehnition of lattice 
topological charge converges to the continuum dehnition when the lattice spacing is sent to 
zero, but it relies on the fact that the gauge helds are smooth. The ultraviolet huctuations 
of the gauge conhgurations are smoothed out via the cooling method nnn before Ql is 
measured. The cooling method is well tested and understood in the context of topological 
charge measurements on the lattice and it is equivalent to a smearing procedure or a gradient 
how smoothing [102] . Empirically it was noted that this smoothing procedure is such that 
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T/Te 

d a^/a Nr No- N^eas 

central value ± statistical error for 

XM XZ Xa Xf 

1.2 

6.001 0.2161 6 64 14000 

0.3880 0.0012 

0.3814 0.0012 

0.3871 0.0012 

0.4192 0.0013 

1.31 

6.053 0.1979 6 48 15600 

0.3495 0.0009 

0.3130 0.0009 

0.3392 0.0010 

0.3691 0.0011 


64 36000 

0.3424 0.0006 

0.3358 0.0006 

0.3402 0.0007 

0.3703 0.0007 


80 14000 

0.3426 0.0010 

0.3389 0.0010 

0.3416 0.0010 

0.3735 0.0011 


6.242 0.1484 8 64 33998 

0.3634 0.0010 

0.3493 0.0010 

0.3520 0.0010 

0.3687 0.0010 


96 14000 

0.3556 0.0015 

0.3533 0.0014 

0.3537 0.0015 

0.3703 0.0015 

1.4 

6.095 0.1852 6 64 54000 

0.3153 0.0005 

0.3077 0.0005 

0.3095 0.0005 

0.3370 0.0005 

1.5 

6.139 0.1729 6 64 54000 

0.2928 0.0005 

0.2833 0.0005 

0.2814 0.0005 

0.3068 0.0005 

1.6 

6.182 0.1621 6 64 53998 

0.2721 0.0005 

0.2587 0.0005 

0.2568 0.0005 

0.2799 0.0005 

1.7 

6.223 0.1525 6 64 24000 

0.2536 0.0008 

0.2330 0.0008 

0.2369 0.0008 

0.2585 0.0008 

1.8 

6.263 0.1441 6 64 24000 

0.2343 0.0008 

0.2005 0.0009 

0.2178 0.0008 

0.2368 0.0008 


80 32000 

0.2320 0.0006 

0.2262 0.0006 

0.2185 0.0006 

0.2368 0.0006 


6.471 0.1080 8 96 14000 

0.2306 0.0016 

0.2170 0.0017 

0.2236 0.0015 

0.2312 0.0016 

1.9 

6.301 0.1365 6 64 24000 

0.2175 0.0009 

0.1672 0.0011 

0.2019 0.0008 

0.2190 0.0009 


80 34000 

0.2164 0.0006 

0.2095 0.0006 

0.2026 0.0006 

0.2189 0.0006 

1.99 

6.550 0.0973 8 64 14795 

0.2013 0.0034 

0.1800 0.0036 

0.1986 0.0029 

0.2013 0.0034 

2.0 

6.338 0.1297 6 48 15600 

0.2040 0.0018 

0.1292 0.0027 

0.1898 0.0016 

0.2042 0.0018 


64 25598 

0.2032 0.0010 

0.1390 0.0014 

0.1893 0.0009 

0.2041 0.0010 


80 26000 

0.2014 0.0008 

0.1920 0.0008 

0.1888 0.0007 

0.2030 0.0008 


96 14000 

0.2004 0.0008 

0.1961 0.0008 

0.1900 0.0008 

0.2038 0.0009 

2.1 

6.373 0.1235 6 80 24000 

0.1880 0.0009 

0.1749 0.0009 

0.1774 0.0008 

0.1889 0.0009 

2.5 

6.502 0.1037 6 128 14000 

0.1497 0.0010 

0.1479 0.0010 

0.1494 0.0008 

0.1492 0.0010 


144 15797 

0.1525 0.0008 

0.1513 0.0008 

0.1495 0.0006 

0.1518 0.0008 


TABLE I: A summary of our lattice parameters and lattice results for the topological susceptibility 
in units of the critical temperature. The temperature in units of the critical temperature and the 
corresponding lattice coupling j3 are related by Eq. (26). We also report the corresponding value 
of the interpolated string tension. The combined number of measurements Nmeas on two different 
streams of configurations, one from a random start and one from a unit start. Our representative 
values at each temperature are emphasized. They are all chosen from the globally-fit definition 
of Q as given in Eq (28d) because it leads to smaller discretization and finite-volume effects. These 
values are going to be used in all following analysis. The different definitions of the topological 
charge Q from which the susceptibility is estimated are described in the text in Eq. (28). 


the multiplicative and an additive renormalization constants to Ql are close to one and 
zero, respectively. This was understood theoretically in Ref. ra. Our use of the cooling 
method should be viewed as a cheaper alternative to the gradient flow method which will 
be employed in our analysis when going beyond the Yang-Mills case. 

Although UV fluctuations are removed by the cooling procedure while keeping topological 
properties largely unchanged, lattice artifacts affect Ql- Since we are interested in continuum 
physics, we use four different dehnitions of the topological charge based on Ql. The following 
dehnitions should all agree in the continuum limit and their behavior as a function of a will 
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Monte Carlo time Frequency 

FIG. 2: The topological charge as a function of Monte Carlo time for the 80^x6 /3 = 6.301 
(T/Tc = 1.90) ensemble. In the left panel, the yellow points are the measurements of Qr and the 
purple points are the globally-fit values Qf. In both panels, the brown dotted lines indicate the 
local maxima for the distribution of Qr. The right panel shows a histogram of those points with 
a bin width of 0.05. The dashed purple line is a gaussian distribution with a standard deviation 
given by the second moment of the Qf distribution. 


help us characterize lattice discretization effects on our hnal results: 


= Ql 

Qz = round(QL) 

Qa = round(QL) — “narrow instantons” 

Qf = round(a;QL) a = min ([q:Ql — round(Q;Q 2 ,)]^) 


(28a) 

(28b) 

(28c) 

(28d) 


where the M, Z, a, and / subscripts indicate the real-valued lattice dehnition of Eq. (27), the 


rounded dehnition, a lattice-artifact corrected dehnition that subtracts contributions from 
narrow instantons nnn, and a globally-ht dehnition that redehnes the charge based on the 
properties of the whole distribution [HSIEH]- We measure Ql every 10 Monte-Carlo updates 
and this gives us an autocorrelation time of 1 measurement or smaller for all W = 6 lattices 
and between 2 and 3 measurements for all W = 8 ensembles. In Fig. |^we compare Qr and 
Qf for T/Tc = 1.90 (/3 = 6.301). 

We dehne the lattice topological susceptibility from each of the dehnitions above: 




= lim 

V^oo 


m 

V 


(29) 


where the index i runs over the dehnitions in Eq. (28a)-(28d) and V = a^{N„)^NT-. We 
dehne our inhnite volume susceptibility by measuring y on diherent spatial volumes and by 
choosing the value on the largest volume that does not show signs of change. In practice 
we hnd that the dehnition Xa has the largest discretization ehects while Xf has negligible 
ehects due to hnite lattice spacing. In the following section we estimate hnite-volume and 
discretization ehects for all value of T/Tc explored. This will allow us to choose a hnal value 
of y at each T/Tc that is free of artifacts and can be htted to models in Sec. 
are summarized in Table HI 


VII All results 
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FIG. 3: A study of finite-volume effects at three different temperatures. All lattices have a 
temporal extent Nj- = 6. Each ensemble is slightly offset from its temperature for ease of visibility. 
The statistical errors are smaller than the markers shown. 
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VI. DISCRETIZATION AND FINITE-VOLUME ERRORS 


In this section we stndy the systematic errors associated with having a discretized hnite 
volume and a hnite lattice spacing. We start by looking at the various dehnitions Xi fixed 
T/Tc and for increasingly larger boxes. When the physical volume is large enough, we expect 
to obtain a constant topological susceptibility corresponding to Eq. (29). For the simple 


rounded definition (shown in red in Fig. |^, the smaller volumes have a consistently 
higher topological susceptibility. In contrast, the artifact-corrected integer definition Xa 
(shown in blue), has consistently smaller topological susceptibilities for smaller volumes. 
The unrounded xm. (yellow) and globally-fit rounded Xf (purple) definitions have essentially 
no finite-volume corrections, but the unrounded measure is systematically lower than the 
globally-fit definition. At high temperature (small lattice spacing), Xf matches xz whereas 
Xm does not. For this reason we choose X/ as our most reliable definition. This comports 
with the findings of Ref. |69]. Because of the concordance between the different volumes it 
seems as though many of the ensembles are all in the infinite-volume limit. It is apparent 
from Fig. I^hat Xf has hardly any finite-volume effects. 

In Fig.^ we show discretization effects at two different temperatures. At T/Tc=1.31 
we show the topological susceptibility for the same physical volume at two different lattice 
spacings. At T/Tc=1.8 we show the topological susceptibility for physical spatial extents 
that differ by roughly 13%. However, since we know from our finite-volume study shown 
in Fig. [^that the volume makes essentially no difference at high temperatures (except for 
the artifact-corrected definition Xa), we can assume the difference in physical volume to 
be negligible compared to the discretization errors. For both temperatures, Xf does not 
change when the lattice spacing is decreased by ~ 30%. For the higher temperature, which 
corresponds to an overall smaller lattice spacing, a similar effect is observed for xz- Because 
the globally-fit definition of the topological susceptibility is essentially independent of the 
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FIG. 4: A study of lattice spacing effects at two different temperatures. At r/Tc=1.31 the lattices 
shown are 48^ x 6 and 64^ x 8 which have exactly the same physical volume, while at T/Tc=1.8 
the lattice shown are 64^ x 6 and 96^ x 8, which have physical spacial extents that differ by only 
13%. Each ensemble is slightly offset from its temperature for ease of visibility. The statistical 
errors are smaller than the markers shown. 

volume and lattice spacing, we use that definition for our final “curated” data set, shown 
in bold in Table In the following we only consider a statistical uncertainty on the data 
points since the systematics effects discussed above are negligible. 
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VII. COMPARISON BETWEEN LATTICE RESULTS AND MODELS 


We £t our best data points, those in bold in Table |T] and shown in Fig. to the forms 


shown in (23). Since the lattice calculations naturally yield x/T^ and our temperatures are 


naturally in units of Tc we actually fit 


Adigm 


c 


{T/T,y 


Aiilm 

2^4 


3^0 


(T/ry 


-di 


exp 


do 


Tr 


+ do, I In — 


(30) 


as functions of T/Tc rather than T itself. This allows us to postpone consideration of the 
systematic error arising from setting the scale in physical units with an uncertain critical 
temperature. 

The central values shown in Fig. and are the results from htting all the data points. 
The £t parameters are shown in Table with a statistical uncertainty of one standard 
deviation (Icr). 

For the systematic fitting error on the DIGM we fit the same form to every subset of the 
best points with cardinality 3 or greater and take the outer envelope of all such hts’ statistical 
Icr error bands. This procedure is extremely conservative. The resulting systematic fitting 
error is shown as a light purple band in Fig. where the statistical error band is entirely 
covered by the width of the line. 
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DIGM 

xVd.o.f. = 1.2 
C 0.0869 ± 0.0015 
n 5.64 ± 0.04 



IILM 


X' 

Vd.o.f 


= 1.7 

gdo 

0.079 

± 

0.006 

di 

-4.9 

± 

0.5 

d2 

-1.7 

± 

1.0 

dz 

1.2 

± 

0.7 


TABLE II: Fit parameters for the DIGM and IILM fit to all of our best data points. 

The DIGM is not reliable at low temperatures, and so that region is not shown in Fig. 
The systematic errors at high temperature decrease—all the hts grow closer to the central 
value. 

In the pure-glue case it is known EH that the DIGM exhibits, at high temperatures, the 
leading behavior y ~ T“"' where 

n = y iV, -4 = 7 (31) 

for SU(3), compared to our numerical result n = 5.64T0.04. This suggests the temperatures 
we studied are not high enough to trust the high-temperature DIGM parameters. 

The IILM is more reliable at low temperatures and unreliable at high temperatures. 
Indeed, our best-£t IILM curve has a turning point at T/Tc ~ 5.8, beyond which it exhibits 
an unphysical increase of the topological susceptibility with temperature. 

Because the IILM is not designed to capture high-temperature dynamics and has 
markedly more free parameters, the enveloping approach we used for estimating the sys¬ 
tematic error bands for the DIGM is far too conservative and gets dominated by hts to 
a few high-temperature points. Instead, we follow a jackknife-inspired procedure: for the 
central value we fit to all the best points, and for the systematic htting error we £t to all 
subsets where one data point is omitted. The outer envelope of all those hts’ Icr statistical 
error bands is shown as the systematic htting error in Fig. The inner, darker error band 
shown there is the statistical error band on the best ht of all of our data points. 


VIII. ILLUSTRATIVE EXAMPLE USING LATTICE DATA TO EXTRACT AX- 
ION MASS BOUNDS 

As discussed in previous sections, the topological susceptibility y for QGD as a function 
of temperature is required as input for the cosmological evolution of axions in the early 
universe. In this section, we will take our lattice results for y and carry out this procedure to 
calculate axion mass bounds with a signihcant disclaimer: our lattice calculation is ahected 
by the absence of fermionic degrees of freedom. This approach means that our numerical 
calculation was significantly cheaper than it would have been if we had calculated with full 
QGD, but can lead to a host of uncontrolled systematics. Most obvious is the diherence 
in critical temperatures between Yang-Mills and QGD^^. The critical temperature is the 


There is roughly a factor of two difference between QCD {T^ ~ 154 MeV [SUIIMIISI]) and Yang-Mills 
(Tc ~ 284 MeV [7H]). 
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FIG. 5: The largest volumes from Gattringer et al. m and our best Nr = 6 data points at each 
temperature. The statistical errors on our points are smaller than the markers shown. 
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FIG. 6: A ht of our best points to the dilute instanton gas model (DIGM). 


primary scale where the theory deconhnes. It also indicates the onset of the fall of x as 
temperature grows. In each of the individual steps to this point, this has not been a issue, 
as all quantities discussed in the result sections have been dimensionless ratios that do not 
depend on scale setting. However, in the hnal steps required to extract a bound, the relative 
size of Tc to the present day temperature of the cosmic microwave background, Ty, plays 
a role. It should be noted that this systematic will not be an issue for future calculations 
with physical QCD parameters. For now, it prevents bounds relevant to reality from being 
extracted. 
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FIG. 7: A fit of our best points to the interacting instanton liquid model (IILM). 


Let us describe the overall effect of setting the absolute scale Tc- As discussed in Sec. |II A 
the axion begins rolling down the potential at the temperature where rua ~ 3iL. With simple 
algebra this may be rewritten 

x(Ti) ^ . (32) 


The fortuitous form of H in Eq. (11), 


H^{T) = 


TT 




g*{T) 


(33) 


makes it particularly nice to write the relation (32) in terms of T/Tc^ the temperature ratio 
that shows up in our lattice calculations. One may show that the relation irta ~ 3i7 can be 
rewritten 

^ (T, • 


A(ri/r.) 

c 


10M2; 


(7.(Ti/Ta-Ta)- ( ^ 


(34) 


which allows us to put all of the absolute scale dependence into the argument of ( 7 *. Fortu¬ 
nately, at the temperatures of interest is insensitive to the difference between the critical 
temperature in QCD and the pure-glue theory 0 — it is essentially a constant at and be¬ 


tween those temperatures. The particular temperature dependence of Eq. (34) also means 


that we do not need to translate our lattice results for into physical units, therefore 
avoiding what is usually the largest source of systematic uncertainties. 

However, this does not imply that all the systematic errors are erased. Instead, because 
( 7 * captures the relevant degrees of freedom of reality, which includes quarks, rather than a 
quark-free universe, we still have an uncontrolled systematic error. This systematic will be 
avoided when future lattice calculations use physical QCD parameters, while the uncertainty 
in the QCD critical temperature Tc will be remain irrelevant thanks to the insensitivity of 
77 *. Thus, we expect the reliability of future lattice calculations to be controlled entirely by 
the accuracy of the determination of x/T^ and cosmological inputs, not the accuracy of the 
QCD critical temperature. 
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FIG. 8; The DIGM-inspired fit with systematic fitting uncertainty shown with lines of {3H)‘^ for 
three representative choices of fa- 


Of course, with x/T^ from QCD in hand, it will make sense to not simply estimate the 
onset of the axion held rolling when nia ~ 3H but to instead solve the cosmological equations 
of motion numerically. Moreover, it will be simple to put reliable bounds on the pre-inhation 
PQ-breaking scenario, constraining a combination of the initial 6 parameter and the axion 
mass. Therefore we emphasize our results in the rest of this section as the hrst step towards 
a full-hedged calculation of axion constraints from hrst principles. 

In Fig. 1^ we plot our numerically-ht DIGM extrapolation of x/T^, the left-hand side 
of Eq. ( [M] ) and three example right-hand sides for diherent choices of /„ as a function of 
temperature. The intersection gives Ti/T^ as a function of fa, which is shown in Fig. In 
terms of the DIGM ht parameters, we have 


Ti 

IOC 

1 

to 

_1 

Tc 


-1 


(35) 


where the insensitivity of to temperatures in the regime of interest allows us to solve 
self-consistently with ease. In a full calculation without extrapolation using a model, Ti can 
be determined as a function of fa numerically from Eq. (34). In Fig. our extracted Ti/Tc 
is compared with that of the IILM [54] . 

With Ti in hand, we can calculate the axion energy density using 


= 


Pa{Ty) 


Pc 


Pc = 3.978 X 10 


-11 




0.701 


eV^ 


(36) 


The resulting energy density from the misalignment mechanism is given by 


h^^a = 0.107(1) 


F{ei)el 

27 


fa 


1012 GeV 


1.2074(8) 


(37) 
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FIG. 9: The value of Ti, the temperature where x = as a function of fa if X is given by 

our DIGM-inspired fits. The purple band includes the systematic fitting uncertainty and the mild 
scale setting ambiguity added in quadrature. The dashed line is the results from Ref. [5l] when 
Ta = 154 MeV. 


where the errors are determined by bootstrapping the fit results for the DIGM shown in 
Table |TT| A bound can be derived for the post-inflation PQ-breaking scenario by starting 
from this density, setting (Of) = vr^/S, and comparing with the latest dark matter density 
from PLANCK [103] 

Qah'^ < ^DMh^ = 0.1199 ± 0.0027. (38) 


This procedure leads to an upper bound on the value of fa and Eq. ([^ can be used to 
translate this into a lower bound on the axion mass today. Putting all these pieces together 
using lattice inputs for Eq. (15), we arrive at the bounds 


Post-inflation PQ breaking, 
physical g^, m^, U, T^, 
pure-glue x 



< (4.10 ±0.04) X 10^1 GeV 
> (14.6 ±0.1) /reV 


(39) 


which includes all the statistical, systematic and scale-setting errors from fitting to and 
extrapolating with the DIGM but does not include any cosmological uncertainties. Also, 
to reiterate, this is only the bound from the misalignment mechanism and does not include 
string or domain wall decay, which would give even tighter bounds. In comparison, the 
most comprehensive, up-to-date bound from the misalignment mechanism is given by rria > 
21 peV [5^ . 

Again, it should be noted that these results are for the topological susceptibility x from a 
pure-glue calculation and not from the full QCD theory with dynamical fermionic degrees of 
freedom. That being said, this calculation is still illuminating and illustrates how uncertainty 
quantification due to strong dynamics can be propagated to the final steps of this calculation. 
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IX. CONCLUSIONS 


We performed a lattice calculation of topological susceptibility for pure Yang-Mills theory 
at high temperatures and applied the results to the problem of axion cosmology. There are 
two primary accomplishments of this work. 

First, we developed the connection between actual lattice data and the well-known pro¬ 
cedure for calculating the axion density, with controlled uncertainties in the input from the 
non-perturbative free energy. Once lattice QCD results span the temperatures of interest, 
a simple interpolating £t can be used to solve the cosmological equations of motion. In the 
mean time, we must rely on model-inspired £t forms extrapolate to the high temperatures 
of interest. In this regard, we find that the DIGM-like £t, Eq. (23), works remarkably well 
for a theory consisting only of gluons, yielding a y^/d.o.f ~ 1.2 with very precise lattice data 
(statistical errors roughly at the 0.2% level). However, we find that our points do not con¬ 
form to the pure-glue expectation that the DIGM exponent n is 7, but that n = 5.64 ±0.04. 
Also, even when propagating the systematic errors, the hnal bounds were found to have total 
uncertainties only at the few percent level. Thus, the DIGM-like £t proves to be a good 
continuous £t form for translating hnite lattice calculations at discrete temperature values 
to an equation that can be used to solve the cosmology evolution equations numerically. 
Moreover, the lattice can provide reliable errors on the £t parameters and these errors can 
be propagated to the hnal answers straightforwardly. 

Second, we aimed to perform a detailed study of lattice systematic effects of lattice 
calculations of the topological susceptibility at high temperature. To accomplish this, small 
statistical errors are a must. Moreover multiple volumes, lattice spacings and discretized 
definitions of the topological susceptibility are required to assess systematic effects and add 
precision. As expected, it was found that the lattice artifacts affect each of the discretized 
definitions of the topological susceptibility differently, and while lattice spacing effects were 
the predominant error in most of the definitions, one definition (artifact-corrected, Qa) was 
extremely sensitive to finite volumes at high temperatures/small lattice spacing, leading to 
shifts of order 100%. Overall, it was found, as in Ref. [HHIES], that the globally-£t dehnition, 
Qj-, had virtually no systematic errors and while a distinct systematic error could be seen 
in all of the other dehnitions, no effect larger than the statistical error bars could be quoted 
for this dehnition. 

One disclaimer that has been pointed out frequently throughout this paper is that the 
lattice results were for a purely gluonic theory and the procedure for determining the axion 
mass bounds assumed full QGD. This leads to multiple issues, which essentially reduce 
the results in this paper to simply be a preliminary step towards being useful in realistic 
experiments. The low temperature axion mass depends on the pion mass and decay constant, 
which do not exist in a purely gluonic theory, and the relevant number of degrees of freedom, 

depend on the number of fermions below the temperatures of interest. Second, QGD 
and pure Yang-Mills diher in deconhnement temperatures by roughly a factor of two when 
scales are hxed to heavy-quark scales (thus, correctly describing high temperature physics 
and heavy-quark physics cannot be done simultaneously). Third, dynamical fermion zero- 
modes are expected to play an important role in the temperature behavior of the topological 
susceptibility, which are omitted in Yang-Mills lattice calculations. 

Despite all of these limitations, it is tempting to see how the bounds behave compared 
to quoted values of the axion mass bounds. The results in Eq. (39) imply that the axion 
mass bound is weaker that the up-to-date value of > 21 peV, and if true, would imply 
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more viable parameter space for the post-inflation PQ-breaking scenario. Current and next 
generation ADMX experiments will throughly explore this additional axion mass parameter 
space over the next few years. Moreover, these lattice-based bounds have a well dehned and 
robust theoretical error associated with them, which is a very compelling feature compared 
to previous bounds. 

While our bounds were derived without the inclusion of fermions in the lattice calcula¬ 
tion, they are ample motivation for future studies and calculations in this direction. From 
this work many of the concerns on volume limitations have been quelled, which suggests 
that lattice QCD calculations on smaller volumes may be sufficient and prohibitively large 
volumes may not be needed. The lattice QCD issue of the increased auto-correlation and 
freezing of topological charge could still be signihcant, but one potential solution is to extract 
the topological susceptibility with hxed topology as in Ref. m or work with anisotropic 
lattices to boost the volume in space while keeping the temporal direction small to go to 
high temperatures. Both of these approaches should be explored further at higher precision 
to accurately quantify the size of the systematic effects. Once methods of extracting topo¬ 
logical susceptibilities at high temperatures are reliable, the next step would be to perform 
this calculation for QCD. Typically, calculations of this kind start with larger-than-physical 
quark masses (i.e. pion masses around 400 MeV) with non-chiral fermion discretizations. 
However, using lattices with physical quark masses and chiral domain-wall fermions [60] is 
most likely next step given current computing resources and given that much of the involved 
scale setting procedure and non-perturbative tuning has been performed. At that stage, 
the lattice QCD results will be directly applicable to cosmological simulations and will have 
direct connection with experimental searches. 
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